home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 19 / CU Amiga Magazine's Super CD-ROM 19 (1998)(EMAP Images)(GB)[!][issue 1998-02].iso / CUCD / Programming / LEDA / source / src / arith / imul.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-11-16  |  11.5 KB  |  582 lines
  1. /*******************************************************************************
  2. +
  3. +  LEDA  3.1c
  4. +
  5. +
  6. +  imul.c
  7. +
  8. +
  9. +  Copyright (c) 1994  by  Max-Planck-Institut fuer Informatik
  10. +  Im Stadtwald, 6600 Saarbruecken, FRG     
  11. +  All rights reserved.
  12. + 
  13. *******************************************************************************/
  14.  
  15.  
  16. /*     imul.c        RD, 29.01.90    */
  17. /*    bug in IasIsrint fixed, RD, 32.08.91    */
  18. /*    sign bug in IasImuP fixed, RD, 26.2.92    */
  19. /*    Karatsuba-Verfahren, RD, 20.7.92    */
  20. /*    bug in memory usage in do_karatsuba, RD, 14.14.92    */
  21.  
  22. #include <LEDA/impl/iint.h>
  23. #include <LEDA/impl/iloc.h>
  24.  
  25.  
  26. #ifdef DO_NOT_USE_KARATSUBA
  27.  
  28. void IasImuI (accu, a, b)
  29.     pInteger accu;
  30.     const Integer *a, *b;
  31. /* accu=a*b; */
  32. {       register int i, al, bl, neededlength;
  33.         register pPLACE paccu, pa, pb;
  34.     if (accu==a) {
  35.         ImuasI(accu, b);
  36.         return;
  37.     }
  38.     if (accu==b) {
  39.         ImuasI(accu, a);
  40.         return;
  41.     }
  42.     al=a->length;
  43.     bl=b->length;
  44.         if (!al || !bl) {
  45.         Iasint(accu, 0);
  46.         return;
  47.         }
  48.         neededlength=al+bl;
  49.     if (neededlength > accu->maxlength) {
  50.         delvec(accu->vec, accu->maxlength);
  51.         accu->maxlength=neededlength;
  52.         paccu=newvec(&accu->maxlength);
  53.         accu->vec=paccu;
  54.     } else
  55.         paccu=accu->vec;
  56.         pa=a->vec;
  57.     pb=b->vec;
  58.     if (al > bl) {
  59.             for ( i = 0; i < al; i++)
  60.                 paccu[i] = 0;
  61.             for ( i = 0; i < bl; i++) 
  62.                 paccu[i+al] = 
  63.             vecmuladd(paccu+i, pa, pb[i], al);
  64.         } else {           
  65.             for ( i = 0; i < bl; i++)
  66.                 paccu[i] = 0;
  67.             for ( i = 0; i < al; i++) 
  68.                 paccu[i+bl] = 
  69.             vecmuladd(paccu+i, pb, pa[i], bl);
  70.         }
  71.         if (paccu[neededlength-1])
  72.             accu->length = neededlength;
  73.         else
  74.             accu->length = neededlength-1;
  75.         accu->sign = a->sign ^ b->sign;
  76. }        /* IasImuI */
  77.  
  78. void ImuasI(a, b)
  79.     pInteger a;
  80.     const Integer *b;
  81. /* a*=b; */
  82. {       int maxl;
  83.         register pPLACE paccu, pa, pb;
  84.         register int i, al, bl;
  85.     al=a->length;
  86.     bl=b->length;
  87.         if (!al || !bl) {
  88.         Iasint(a, 0);
  89.         return;
  90.         }
  91.         maxl=al+bl;
  92.         paccu=newvec(&maxl);
  93.         pa=a->vec;
  94.     pb=b->vec;
  95.     if (al > bl) {
  96.             for ( i = 0; i < al; i++)
  97.                 paccu[i] = 0;
  98.             for ( i = 0; i < bl; i++) 
  99.                 paccu[i+al] = 
  100.             vecmuladd(paccu+i, pa, pb[i], al);
  101.         } else {           
  102.             for ( i = 0; i < bl; i++)
  103.                 paccu[i] = 0;
  104.             for ( i = 0; i < al; i++) 
  105.                 paccu[i+bl] = 
  106.             vecmuladd(paccu+i, pb, pa[i], bl);
  107.         }
  108.     al+=bl;
  109.         if (paccu[al-1])
  110.             a->length = al;
  111.         else
  112.             a->length = al-1;
  113.         a->sign ^= b->sign;
  114.     delvec(a->vec, a->maxlength);
  115.     a->vec=paccu;
  116.     a->maxlength=maxl;
  117. }        /* ImuasI */
  118.  
  119. #else
  120. /* use karatsuba */
  121.  
  122. /* Prerequisite: enough space in the vectors to increase the length
  123.     appropriately, given with the built-in memory management.
  124. */
  125.  
  126. #ifndef KARATSUBA_LIMIT1
  127. #define KARATSUBA_LIMIT1 15
  128. #endif
  129. #ifndef KARATSUBA_LIMIT2
  130. #define KARATSUBA_LIMIT2 KARATSUBA_LIMIT1
  131. #endif
  132.  
  133. static void do_karatsuba( /* paccu, pa, pb, al, bl */ );
  134.  
  135. void IasImuI (accu, a, b)
  136.     Integer *accu;
  137.     const Integer *a, *b;
  138. /* accu=a*b; */
  139. /* Karatsuba, verwendet do_karatsuba */
  140. {       register int i, al, bl, neededlength;
  141.         register pPLACE paccu, pa, pb;
  142.  
  143.     if (accu==a) {
  144.         ImuasI(accu, b);
  145.         return;
  146.     }
  147.     if (accu==b) {
  148.         ImuasI(accu, a);
  149.         return;
  150.     }
  151.     al=a->length;
  152.     bl=b->length;
  153.         if (!al || !bl) {
  154.         Iasint(accu, 0);
  155.         return;
  156.         }
  157.         neededlength=al+bl;
  158.     if (neededlength > accu->maxlength) {
  159.         delvec(accu->vec, accu->maxlength);
  160.         accu->maxlength=neededlength;
  161.         paccu=newvec(&accu->maxlength);
  162.         accu->vec=paccu;
  163.     } else
  164.         paccu=accu->vec;
  165.     if (al > bl) {
  166.         if (bl <= KARATSUBA_LIMIT1) {
  167.             /* Standardmultiplikation */
  168.             pa=a->vec;
  169.         pb=b->vec;
  170.         for ( i = 0; i < al; i++)
  171.             paccu[i] = 0;
  172.         for ( i = 0; i < bl; i++) 
  173.             paccu[i+al] = 
  174.                 vecmuladd(paccu+i, pa, pb[i], al);
  175.         } else
  176.         do_karatsuba(paccu, a->vec, b->vec, al, bl);
  177.         } else {
  178.         if (al <= KARATSUBA_LIMIT1) {       
  179.             /* Standardmultiplikation */
  180.             pa=a->vec;
  181.         pb=b->vec;
  182.         for ( i = 0; i < bl; i++)
  183.             paccu[i] = 0;
  184.         for ( i = 0; i < al; i++) 
  185.             paccu[i+bl] = 
  186.                 vecmuladd(paccu+i, pb, pa[i], bl);
  187.         } else
  188.         do_karatsuba(paccu, b->vec, a->vec, bl, al);
  189.         }
  190.         if (paccu[neededlength-1])
  191.             accu->length = neededlength;
  192.         else
  193.             accu->length = neededlength-1;
  194.         accu->sign = a->sign ^ b->sign;
  195. }    /* IasImuI, karatsuba    */
  196.  
  197. void ImuasI(a, b)
  198.     pInteger a;
  199.     const Integer *b;
  200. /* a*=b; */
  201. /* Karatsuba, verwendet do_karatsuba */
  202. {       int maxl;
  203.         register pPLACE paccu, pa, pb;
  204.         register int i, al, bl;
  205.     al=a->length;
  206.     bl=b->length;
  207.         if (!al || !bl) {
  208.         Iasint(a, 0);
  209.         return;
  210.         }
  211.         maxl=al+bl;
  212.         paccu=newvec(&maxl);
  213.     if (al > bl) {
  214.         if (bl <= KARATSUBA_LIMIT1) {
  215.             /* Standardmultiplikation */
  216.             pa=a->vec;
  217.         pb=b->vec;
  218.         for ( i = 0; i < al; i++)
  219.             paccu[i] = 0;
  220.         for ( i = 0; i < bl; i++) 
  221.             paccu[i+al] = 
  222.                 vecmuladd(paccu+i, pa, pb[i], al);
  223.         } else
  224.         do_karatsuba(paccu, a->vec, b->vec, al, bl);
  225.         } else {
  226.         if (al <= KARATSUBA_LIMIT1) {       
  227.             /* Standardmultiplikation */
  228.             pa=a->vec;
  229.         pb=b->vec;
  230.         for ( i = 0; i < bl; i++)
  231.             paccu[i] = 0;
  232.         for ( i = 0; i < al; i++) 
  233.             paccu[i+bl] = 
  234.                 vecmuladd(paccu+i, pb, pa[i], bl);
  235.         } else
  236.         do_karatsuba(paccu, b->vec, a->vec, bl, al);
  237.         }
  238.     al+=bl;
  239.         if (paccu[al-1])
  240.             a->length = al;
  241.         else
  242.             a->length = al-1;
  243.         a->sign ^= b->sign;
  244.     delvec(a->vec, a->maxlength);
  245.     a->vec=paccu;
  246.     a->maxlength=maxl;
  247. }        /* ImuasI, karatsuba */
  248.  
  249.  
  250. static void do_karatsuba(paccu, pa, pb, al, bl)
  251.     PLACE *paccu, *pa, *pb;
  252.     int    al, bl;
  253.     /* now pb is the shorter number */
  254. {        PLACE *tprod, *tmp, *res;
  255.     int i, tmpl, k, q, r, n, n0, resl;
  256.  
  257.     /* calculate k minimal with 
  258.         bl <= 2^k*n0 and n0<=KARATSUBA_LIMIT2, n=2^k*n0 */
  259.     n0=bl;
  260.     k=0;
  261.     while (n0>KARATSUBA_LIMIT2) {
  262.         n0=(n0+1)/2;
  263.         k++;
  264.     }
  265.     n=n0<<k;
  266.     /* fill a, b with zeroes, maybe al <= n */
  267.     for (i=bl; i<n; i++)
  268.         pb[i]=0;
  269.     for (i=al; i<n; i++)
  270.         pa[i]=0;
  271.     /* split a */
  272.     q=al/n;
  273.     r=al%n;
  274.     if (q==0) {    /* case al < n */
  275.         q=1;
  276.         r=0;
  277.     }
  278.     /* Now get temporaries for result of karatsuba_mult and tmp */
  279.     tmpl=2*n;
  280.     tmp=newvec(&tmpl);
  281.     tprod=newvec(&tmpl);
  282.     /* Get temporary for result. */
  283.     resl=al+n;
  284.     res=newvec(&resl);
  285.     /* calculate via karatsuba_mult a_i * b, 0<= i <q */
  286.     karatsuba_mult(res, pa, pb, tmp, n);
  287.     for (i=1; i<q; i++) {
  288.         karatsuba_mult(tprod, &pa[i*n], pb, tmp, n);
  289.         cvadd(&res[i*n], &res[i*n], tprod, n , 2*n);
  290.     }
  291.     /* calculate a_q * b */
  292.     for (i=q*n; i<al; i++)
  293.         res[i+bl] = vecmuladd(&res[i], pb, pa[i], bl);
  294.     for (i=0; i<al+bl; i++)
  295.         paccu[i]=res[i];
  296.     delvec(tmp, tmpl);
  297.     delvec(tprod, tmpl);
  298.     delvec(res, resl);
  299. }        /* do_karatsuba */
  300.  
  301. #endif
  302. /* DO_NOT_USE_KARATSUBA */
  303.  
  304. #ifdef __GNUC__
  305. void IasImuP(Integer *accu, const Integer *b, PLACE c)
  306. #else
  307. void IasImuP(accu, b, c)
  308.     register Integer *accu;
  309.     const Integer *b;
  310.     PLACE c;
  311. #endif
  312. {    register int nl;
  313.     register PLACE carry;
  314.  
  315.     if (accu==b) {
  316.         ImuasP(accu, c);
  317.         return;
  318.     }
  319.     if (!c) {
  320.         Iasint(accu, 0);
  321.         return;
  322.     }
  323.     nl=b->length+1;
  324.     if (nl > accu->maxlength) {
  325.         delvec(accu->vec, accu->maxlength);
  326.         accu->maxlength=nl;
  327.         accu->vec=newvec(&accu->maxlength);
  328.     }
  329.     carry=vecmul(accu->vec, b->vec, c, b->length);
  330.     accu->vec[nl-1]=carry;
  331.     if (carry)
  332.         accu->length = nl;
  333.     else
  334.         accu->length = nl-1;
  335.     accu->sign=b->sign;
  336. }        /* IasImuP */
  337.  
  338. #ifdef __GNUC__
  339. void ImuasP(Integer *accu, PLACE b)
  340. #else
  341. void ImuasP(accu, b)
  342.     register pInteger accu;
  343.     PLACE b;
  344. #endif
  345. {    register PLACE carry, *paccu;
  346.     register int nl;
  347.     BOOLEAN neednew;
  348.     int maxl;
  349.  
  350.     if (!b) {
  351.         Iasint(accu, 0);
  352.         return;
  353.     }
  354.     nl=accu->length+1;
  355.     if (nl>accu->maxlength) {
  356.         neednew=TRUE;
  357.         maxl=nl;
  358.         paccu=newvec(&maxl);
  359.     } else {
  360.         paccu=accu->vec;
  361.         neednew=FALSE;
  362.     }
  363.     carry=vecmul(paccu, accu->vec, b, accu->length);
  364.     paccu[nl-1]=carry;
  365.     if (carry)
  366.         accu->length=nl;
  367.     if (neednew) {
  368.         delvec(accu->vec, accu->maxlength);
  369.         accu->vec=paccu;
  370.         accu->maxlength=maxl;
  371. }    }    /* ImuasP */
  372.  
  373. /****************************************************/
  374.  
  375. void IasIsrint(a, b, count)
  376.     register pInteger a;
  377.     register const Integer *b;
  378.     unsigned int count;
  379. /* Shift nach rechts */
  380. {    register PLACE accu, help, *pa, *pb;
  381.     register int i;
  382.     int pts, bts, bleft, nl;
  383.     if (a==b) {
  384.         Israsint(a, count);
  385.         return;
  386.     }
  387.     pts=count/LOGPLACE;
  388.     if (pts>=b->length) {
  389.         Ias0(a);
  390.         return;
  391.     }
  392.     bts=count%LOGPLACE;
  393.     bleft=LOGPLACE-bts;
  394.     nl=b->length-pts;
  395.     if (nl>a->maxlength) {
  396.         delvec(a->vec, a->maxlength);
  397.         a->maxlength=nl;
  398.         a->vec=newvec(&a->maxlength);
  399.     }
  400.     pa=a->vec;
  401.     pb=b->vec;
  402.     if ( !bts ) {
  403.         for (i=pts; i<b->length; i++)
  404.             pa[i-pts]=pb[i];
  405.         a->length=nl;
  406.         a->sign=b->sign;
  407.         return;
  408.     }
  409.     accu=pb[pts];
  410.     accu>>=bts;
  411.     for (i=pts+1; i<b->length; i++) {
  412.         help=pb[i];
  413.         accu|=(help<<bleft);
  414.         pa[i-pts-1]=accu;
  415.         accu=help>>bts;
  416.     }
  417.     pa[nl-1]=accu;
  418.     if (accu)
  419.         a->length=nl;
  420.     else
  421.         a->length=nl-1;
  422.     if (a->length)
  423.         a->sign=b->sign;
  424.     else
  425.         a->sign=PLUS;
  426. }        /* IasIsrint */
  427.  
  428. void Israsint(a, count)
  429.     register pInteger a;
  430.     unsigned int count;
  431. /* Shift nach rechts */
  432. {    register PLACE accu, help, *p;
  433.     register int i;
  434.     int pts, bts, bleft, l;
  435.     pts=count/LOGPLACE;
  436.     if (pts>=a->length) {
  437.         Ias0(a);
  438.         return;
  439.     }
  440.     bts=count%LOGPLACE;
  441.     bleft=LOGPLACE-bts;
  442.     p=a->vec;
  443.     if ( !bts ) {
  444.         for (i=pts; i<a->length; i++)
  445.             p[i-pts]=p[i];
  446.         a->length=a->length-pts;
  447.         return;
  448.     }
  449.     accu=p[pts];
  450.     accu>>=bts;
  451.     for (i=pts+1; i<a->length; i++) {
  452.         help=p[i];
  453.         accu|=(help<<bleft);
  454.         p[i-pts-1]=accu;
  455.         accu=help>>bts;
  456.     }
  457.     l=a->length-pts;
  458.     if (accu) {
  459.         p[l-1]=accu;
  460.         a->length=l;
  461.     } else
  462.         a->length=l-1;
  463.     if (!a->length)
  464.         a->sign=PLUS;
  465. }        /* Israsint */
  466.  
  467. void IasIslint(a, b, count)
  468.     register pInteger a;
  469.     register const Integer *b;
  470.     unsigned int count;
  471. /* Shift nach links */
  472. {    register PLACE accu, help, *pa, *pb;
  473.     register int i;
  474.     int pts, bts, bleft, nl;
  475.     if (!b->length) {
  476.         Iasint(a, 0);
  477.         return;
  478.     }
  479.     if (a==b) {
  480.         Islasint(a, count);
  481.         return;
  482.     }
  483.     pts=count/LOGPLACE;
  484.     bts=count%LOGPLACE;
  485.     bleft=LOGPLACE-bts;
  486.     nl=b->length+pts+1;
  487.     if (nl>a->maxlength) {
  488.         delvec(a->vec, a->maxlength);
  489.         a->maxlength=nl;
  490.         a->vec=newvec(&a->maxlength);
  491.     }
  492.     a->sign=b->sign;
  493.     pa=a->vec;
  494.     pb=b->vec;
  495.     for (i=0; i<pts; i++) 
  496.         pa[i]=0;
  497.     if ( !bts ) {
  498.         for (i=0; i<b->length; i++)
  499.             pa[i+pts]=pb[i];
  500.         a->length=nl-1;
  501.         return;
  502.     }
  503.     accu=0;
  504.     for (i=0; i<b->length; i++) {
  505.         help=pb[i];
  506.         accu|=(help<<bts);
  507.         pa[i+pts]=accu;
  508.         accu=help>>bleft;
  509.     }
  510.     pa[nl-1]=accu;
  511.     if (accu)
  512.         a->length=nl;
  513.     else
  514.         a->length=nl-1;
  515. }        /* IasIslint */
  516.  
  517. void Islasint(a, count)
  518.     register pInteger a;
  519.     unsigned int count;
  520. /* Shift nach links */
  521. {    register PLACE accu, help, *pa, *pb;
  522.     register int i;
  523.     int pts, bts, bleft, nl, maxl;
  524.     if (!a->length)
  525.         return;
  526.     pts=count/LOGPLACE;
  527.     nl=a->length+pts+1;
  528.     if (nl>a->maxlength) {
  529.         pb=a->vec;
  530.         maxl=nl;
  531.         pa=newvec(&maxl);
  532.     } else {
  533.         pa=pb=a->vec;
  534.     }
  535.     bts=count%LOGPLACE;
  536.     bleft=LOGPLACE-bts;
  537.     accu=0;
  538.     if ( bts) {
  539.         for (i=a->length-1; i>=0; i--) {
  540.             help=pb[i];
  541.             accu|=(help>>bleft);
  542.             pa[i+pts+1]=accu;
  543.             accu=help<<bts;
  544.         }
  545.         pa[pts]=accu;
  546.     } else {
  547.         pa[nl-1]=0;
  548.         for (i=a->length-1; i>=0; i--)
  549.             pa[i+pts]=pb[i];
  550.     }
  551.     for (i=pts-1; i>=0; i--)
  552.         pa[i]=0;
  553.     if (nl>a->maxlength) {
  554.         delvec(a->vec, a->maxlength);
  555.         a->vec=pa;
  556.         a->maxlength=maxl;
  557.     }
  558.     if (pa[nl-1])
  559.         a->length=nl;
  560.     else
  561.         a->length=nl-1;
  562. }        /* Islasint */
  563.  
  564. BOOLEAN Isr1(a)
  565.     register pInteger a;
  566. {    register int l;
  567.     register BOOLEAN b;
  568.     l=a->length;
  569.     if (!l)
  570.         return 0;
  571.     b=vecsr1(a->vec, l);
  572.     l--;
  573.     if (!a->vec[l])
  574.         a->length=l;
  575.     return b;
  576. }        /* Isr1 */
  577.  
  578. BOOLEAN Ieven(a)
  579.     register const Integer *a;
  580. {    return (!((*(a->vec))&1));
  581. }    /* Ieven */
  582.